### Set working directory ### 
setwd("C:/Users/shkim/Dropbox/Brokerage/Replication") 


### Import some R packages ### 
library(data.table) 
library(lubridate) 
library(stringr) 
library(igraph)  
library(lfe)
library(plm)
library(lmtest)
library(sandwich)
library(stargazer)
# trace(stargazer:::.stargazer.wrap, edit = TRUE)


################ 
### Figure 1 ### 
################

### Import data ### 
rm(list = ls())
trades <- fread("pseudo_data_block_trades.csv")

### Panel A ###
plot(trades$date, trades$p_ratio_p01*100, type = "l", pch=19, col="blue", xlab="Date", ylab="", 
     main = "Block volume, as % of off-exchange volume", ylim = c(0, 50), lty = 2, lwd = 2)

lines(trades$date, trades$p_ratio_p05*100, pch=18, col="red", type="l", lwd = 2)

legend("topright", legend=c("1% of ADV", "5% of ADV"), col=c("blue", "red"), lty=c(2,1), lwd = 2)

### Panel B ###
plot(trades$date, trades$p_vol_p01*100, type = "l", pch=19, col="blue", xlab="Date", ylab="", 
     main = "Block volume, as % of CRSP volume", ylim = c(0, 25), lty = 2, lwd = 2)

lines(trades$date, trades$p_vol_p05*100, pch=18, col="red", type="l", lwd = 2)

legend("topright", legend=c("1% of ADV", "5% of ADV"), col=c("blue", "red"), lty=c(2,1), lwd = 2)

### Panel C ### 
options(scipen = 999)

plot(trades$date, trades$avg_block_size_p01/1e+3, type = "l", pch=19, col="blue", xlab="Date", ylab="", 
     main = "Average block size (in thousand shares)", ylim = c(0, 120), lty = 2, lwd = 2)

lines(trades$date, trades$avg_block_size_p05/1e+3, pch=18, col="red", type="l", lwd = 2)

legend("topleft", legend=c("1% of ADV", "5% of ADV"), col=c("blue", "red"), lty=c(2,1), lwd = 2)

### Panel D ### 
options(scipen = 999)

plot(trades$date, trades$avg_block_dvol_p01/1e+6, type = "l", pch=19, col="blue", xlab="Date", ylab="", 
     main = "Average block size (in million dollars)", ylim = c(0, 7), lty = 2, lwd = 2)

lines(trades$date, trades$avg_block_dvol_p05/1e+6, pch=18, col="red", type="l", lwd = 2)

legend("topleft", legend=c("1% of ADV", "5% of ADV"), col=c("blue", "red"), lty=c(2,1), lwd = 2)


################
### Figure 2 ### 
################

### Import data ### 
rm(list = ls())
connections <- fread("pseudo_data_networks.csv")

g <- graph_from_data_frame(connections[, .(mgmt_cd, broker_cd, b_connect)], directed = FALSE)

V(g)$type <- nchar(V(g)$name) == 5
E(g)$weight <- E(g)$b_connect / 100

V(g)$size[!V(g)$type] <- merge(unique(connections[, .(mgmt_cd, m_cent)]), data.table(mgmt_cd = as.integer(V(g)$name[!V(g)$type]), m_order = 1:uniqueN(connections$mgmt_cd)), by = "mgmt_cd")[order(m_order)]$m_cent
V(g)$size[V(g)$type] <- merge(unique(connections[, .(broker_cd, b_cent)]), data.table(broker_cd = as.integer(V(g)$name[V(g)$type]), b_order = 1:uniqueN(connections$broker_cd)), by = "broker_cd")[order(b_order)]$b_cent

V(g)$color[!V(g)$type] <- "mediumblue"
V(g)$color[V(g)$type] <- "red"

V(g)$shape <- "sphere"
V(g)$shape[V(g)$type] <- "sphere"

V(g)$label<- NA

colrs <- c("mediumblue", "red")

V(g)$size[!V(g)$type] <- merge(unique(connections[, .(mgmt_cd, m_cent)]), data.table(mgmt_cd = as.integer(V(g)$name[!V(g)$type]), m_order = 1:uniqueN(connections$mgmt_cd)), by = "mgmt_cd")[order(m_order)]$m_cent
V(g)$size[V(g)$type] <- 1.2 * merge(unique(connections[, .(broker_cd, b_cent)]), data.table(broker_cd = as.integer(V(g)$name[V(g)$type]), b_order = 1:uniqueN(connections$broker_cd)), by = "broker_cd")[order(b_order)]$b_cent
V(g)$size <- 2 + V(g)$size * 6

plot.igraph(g, layout=layout.fruchterman.reingold)


########################
### Table 1, Panel A ###
########################

### Import data ### 
rm(list = ls())
managers <- fread("pseudo_data_summary_mutual_funds_01.csv")

### Report log(aggregate commissions) ### 
res1 <- managers[, .(var_avg = mean(log(agg_comms), na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report number of brokers ### 
rm(res1, res2, fit1, fit2)

res1 <- managers[, .(var_avg = mean(n_brokers, na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report broker's centrality  ### 
rm(res1, res2, fit1, fit2)

res1 <- managers[, .(var_avg = mean(broker_cent, na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


### Import data ### 
rm(list = ls())
funds <- fread("pseudo_data_summary_mutual_funds_02.csv")

### Report log(TNA) ### 
res1 <- funds[, .(var_avg = mean(log(tna), na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report expense ratio ### 
rm(res1, res2, fit1, fit2)

res1 <- funds[, .(var_avg = mean(exp_ratio, na.rm = TRUE) * 100), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report turnover ratio ### 
rm(res1, res2, fit1, fit2)

res1 <- funds[, .(var_avg = mean(turn_ratio, na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


### Import data ### 
rm(list = ls())
holdings <- fread("pseudo_data_summary_mutual_funds_03.csv")

### Report log(market cap) ### 
res1 <- holdings[, .(var_avg = mean(logmktcap, na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report book-to-market ### 
rm(res1, res2, fit1, fit2)

res1 <- holdings[, .(var_avg = mean(btm, na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report momentum ### 
rm(res1, res2, fit1, fit2)

res1 <- holdings[, .(var_avg = mean(mom, na.rm = TRUE)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


########################
### Table 1, Panel B ###
########################

### Read the data ### 
rm(list = ls())
trades <- fread("pseudo_data_summary_abel_noser_01.csv")

### Report number of brokers ###
res1 <- trades[, .(var_avg = mean(log(agg_comms), na.rm = TRUE)), by = .(m_cent_q, date)][order(m_cent_q, date)]

res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report number of brokers ###
rm(res1, res2, fit1, fit2)

res1 <- trades[, .(var_avg = mean(n_brokers, na.rm = TRUE)), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report broker centrality ###
rm(res1, res2, fit1, fit2)

res1 <- trades[, .(var_avg = mean(broker_cent, na.rm = TRUE)), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


### Import data ### 
rm(list = ls())
trades <- fread("pseudo_data_summary_abel_noser_02.csv")

### Report order size / average daily volume ###
res1 <- trades[, .(var_avg = mean(scaled_ticket_volume, na.rm = TRUE)), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report log(average daily volume) ###
rm(res1, res2, fit1, fit2)

res1 <- trades[, .(var_avg = mean(log_avg_volume, na.rm = TRUE)), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report inverse price ###
rm(res1, res2, fit1, fit2)

res1 <- trades[, .(var_avg = mean(inv_prc, na.rm = TRUE)), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report NYSE dummy ###
rm(res1, res2, fit1, fit2)

res1 <- trades[, .(var_avg = mean(nyse, na.rm = TRUE)), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report stock volatility ###
rm(res1, res2, fit1, fit2)

res1 <- trades[, .(var_avg = mean(stock_volatility, na.rm = TRUE) * 100), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report prev. day's return ###
rm(res1, res2, fit1, fit2)

res1 <- trades[, .(var_avg = mean(ret_lag, na.rm = TRUE) * 100), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "var_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


###############
### Table 2 ###
###############

### Read the data ### 
rm(list = ls())
returns <- fread("pseudo_data_retgap.csv")

### Report return gap ### 
res1 <- returns[, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


############### 
### Table 3 ###
############### 

### Import data ### 
rm(list = ls())
trades <- fread("pseudo_data_shortfall.csv")

### Report trading costs ### 
res1 <- trades[, .(es_avg = mean(abn_trade_cost, na.rm = TRUE) * 100), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "es_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


###############
### Table 4 ###
###############

### Import data ### 
rm(list = ls())
returns <- fread("pseudo_data_retgap_flow.csv")

### Report return gap for fund flow Q1 ### 
res1 <- returns[flow_q == 1, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100, n_funds = uniqueN(wficn)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report return gap for fund flow Q2 ### 
rm(res1, res2, fit1, fit2)

res1 <- returns[flow_q == 2, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100, n_funds = uniqueN(wficn)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report return gap for fund flow Q3 ### 
rm(res1, res2, fit1, fit2)

res1 <- returns[flow_q == 3, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100, n_funds = uniqueN(wficn)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report return gap for fund flow Q4 ### 
rm(res1, res2, fit1, fit2)

res1 <- returns[flow_q == 4, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100, n_funds = uniqueN(wficn)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report return gap for fund flow Q5 ### 
rm(res1, res2, fit1, fit2)

res1 <- returns[flow_q == 5, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100, n_funds = uniqueN(wficn)), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


###############
### Table 5 ###
###############

### Import data ### 
rm(list = ls())
returns <- fread("pseudo_data_retgap_family.csv")

### Report return gap for family size Q1 ### 
res1 <- returns[fam_size_q == 1, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report return gap for family size Q2 ### 
rm(res1, res2, fit1, fit2)

res1 <- returns[fam_size_q == 2, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


###############
### Table 6 ###
###############

### Import data ### 
rm(list = ls())
returns <- fread("pseudo_data_retgap_illiquidity.csv")

### Report return gap for holdings illiquidity Q1 ### 
res1 <- returns[illiq_q == 1, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Report return gap for holdings illiquidity Q2 ### 
rm(res1, res2, fit1, fit2)

res1 <- returns[illiq_q == 2, .(retgap_avg = mean(retgap, na.rm = TRUE) * 100), by = .(cent_q, date)][order(cent_q, date)]
res2 <- dcast(res1, date ~  cent_q, value.var = "retgap_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


###############
### Table 7 ###
###############

### Import data ### 
rm(list = ls())
trades <- fread("pseudo_data_shortfall_small_orders.csv")

### Report trading costs for small orders ### 
res1 <- trades[, .(es_avg = mean(abn_trade_cost, na.rm = TRUE) * 100), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "es_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Import data ### 
rm(list = ls())
trades <- fread("pseudo_data_shortfall_large_orders.csv")

### Report trading costs for large orders ### 
res1 <- trades[, .(es_avg = mean(abn_trade_cost, na.rm = TRUE) * 100), by = .(m_cent_q, date)][order(m_cent_q, date)]
res2 <- dcast(res1, date ~  m_cent_q, value.var = "es_avg")
names(res2) <- c("date", paste0("Q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(Q1 ~ 1, data = res2)
fit1[[2]] <- lm(Q2 ~ 1, data = res2)
fit1[[3]] <- lm(Q3 ~ 1, data = res2)
fit1[[4]] <- lm(Q4 ~ 1, data = res2)
fit1[[5]] <- lm(Q5 ~ 1, data = res2)
fit1[[6]] <- lm((Q5 - Q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


############### 
### Table 8 ### 
############### 

### Import data ### 
rm(list = ls())
trades <- fread("pseudo_data_shortfall_dual.csv")

### Report trading costs ### 
res1 <- trades[, .(es_avg = mean(abn_trade_cost, na.rm = TRUE) * 100), by = .(b_cent_q, m_cent_q, date)][order(b_cent_q, m_cent_q, date)]
res1[, b_cent_q := paste0("b", b_cent_q)]
res1[, m_cent_q := paste0("m", m_cent_q)]

res10 <- dcast(res1, date ~ b_cent_q + m_cent_q, value.var = "es_avg")

### Broker quintile 1 ### 
fit11 <- NULL
fit11[[1]] <- lm(b1_m1 ~ 1, data = res10)
fit11[[2]] <- lm(b1_m2 ~ 1, data = res10)
fit11[[3]] <- lm(b1_m3 ~ 1, data = res10)
fit11[[4]] <- lm(b1_m4 ~ 1, data = res10)
fit11[[5]] <- lm(b1_m5 ~ 1, data = res10)
fit11[[6]] <- lm((b1_m5 - b1_m1) ~ 1, data = res10)

fit12 <- NULL
for (i in 1:6) {
  fit12[[i]] <- coeftest(fit11[[i]], vcov.=NeweyWest(fit11[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit12, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Broker quintile 2 ### 
fit13 <- NULL
fit13[[1]] <- lm(b2_m1 ~ 1, data = res10)
fit13[[2]] <- lm(b2_m2 ~ 1, data = res10)
fit13[[3]] <- lm(b2_m3 ~ 1, data = res10)
fit13[[4]] <- lm(b2_m4 ~ 1, data = res10)
fit13[[5]] <- lm(b2_m5 ~ 1, data = res10)
fit13[[6]] <- lm((b2_m5 - b2_m1) ~ 1, data = res10)

fit14 <- NULL
for (i in 1:6) {
  fit14[[i]] <- coeftest(fit13[[i]], vcov.=NeweyWest(fit13[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit14, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Broker quintile 3 ### 
fit15 <- NULL
fit15[[1]] <- lm(b3_m1 ~ 1, data = res10)
fit15[[2]] <- lm(b3_m2 ~ 1, data = res10)
fit15[[3]] <- lm(b3_m3 ~ 1, data = res10)
fit15[[4]] <- lm(b3_m4 ~ 1, data = res10)
fit15[[5]] <- lm(b3_m5 ~ 1, data = res10)
fit15[[6]] <- lm((b3_m5 - b3_m1) ~ 1, data = res10)

fit16 <- NULL
for (i in 1:6) {
  fit16[[i]] <- coeftest(fit15[[i]], vcov.=NeweyWest(fit15[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit16, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Broker quintile 4 ### 
fit17 <- NULL
fit17[[1]] <- lm(b4_m1 ~ 1, data = res10)
fit17[[2]] <- lm(b4_m2 ~ 1, data = res10)
fit17[[3]] <- lm(b4_m3 ~ 1, data = res10)
fit17[[4]] <- lm(b4_m4 ~ 1, data = res10)
fit17[[5]] <- lm(b4_m5 ~ 1, data = res10)
fit17[[6]] <- lm((b4_m5 - b4_m1) ~ 1, data = res10)

fit18 <- NULL
for (i in 1:6) {
  fit18[[i]] <- coeftest(fit17[[i]], vcov.=NeweyWest(fit17[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit18, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Broker quintile 5 ### 
fit19 <- NULL
fit19[[1]] <- lm(b5_m1 ~ 1, data = res10)
fit19[[2]] <- lm(b5_m2 ~ 1, data = res10)
fit19[[3]] <- lm(b5_m3 ~ 1, data = res10)
fit19[[4]] <- lm(b5_m4 ~ 1, data = res10)
fit19[[5]] <- lm(b5_m5 ~ 1, data = res10)
fit19[[6]] <- lm((b5_m5 - b5_m1) ~ 1, data = res10)

fit20 <- NULL
for (i in 1:6) {
  fit20[[i]] <- coeftest(fit19[[i]], vcov.=NeweyWest(fit19[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit20, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

### Broker quintile 5 minus quintile 1 ### 
fit21 <- NULL
fit21[[1]] <- lm((b5_m1 - b1_m1) ~ 1, data = res10)
fit21[[2]] <- lm((b5_m2 - b1_m2) ~ 1, data = res10)
fit21[[3]] <- lm((b5_m3 - b1_m3) ~ 1, data = res10)
fit21[[4]] <- lm((b5_m4 - b1_m4) ~ 1, data = res10)
fit21[[5]] <- lm((b5_m5 - b1_m5) ~ 1, data = res10)

fit22 <- NULL
for (i in 1:5) {
  fit22[[i]] <- coeftest(fit21[[i]], vcov.=NeweyWest(fit21[[i]], lag=5, adjust=TRUE, verbose=TRUE, prewhite = FALSE))
}

stargazer(fit22, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)


###############
### Table 9 ###
###############

### Panel A ### 
rm(list = ls())
funds <- fread("pseudo_data_lehman_centrality.csv")

fit1 <- NULL
fit1[[1]] <- felm(std_cent ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + lehman_q * post | 0 | 0 | mgmt_cd, data = funds)
fit1[[2]] <- felm(std_cent ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + lehman_q * post | date | 0 | mgmt_cd, data = funds)
fit1[[3]] <- felm(std_cent ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + I(lehman_q == 5) * post + I(lehman_q == 4) * post + I(lehman_q == 3) * post + I(lehman_q == 2) * post | 0 | 0 | mgmt_cd, data = funds)
fit1[[4]] <- felm(std_cent ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + I(lehman_q == 5) * post + I(lehman_q == 4) * post + I(lehman_q == 3) * post + I(lehman_q == 2) * post | date | 0 | mgmt_cd, data = funds)

stargazer(fit1, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE, keep.stat = c("n", "adj.rsq"))

### Panel B ### 
rm(list = ls())
returns <- fread("pseudo_data_lehman_retgap.csv")

fit2 <- NULL
fit2[[1]] <- felm(I(retgap*100) ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + lehman_q * post | 0 | 0 | mgmt_cd, data = returns)
fit2[[2]] <- felm(I(retgap*100) ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + lehman_q * post | date | 0 | mgmt_cd, data = returns)
fit2[[3]] <- felm(I(retgap*100) ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + I(lehman_q == 5) * post + I(lehman_q == 4) * post + I(lehman_q == 3) * post + I(lehman_q == 2) * post | 0 | 0 | mgmt_cd, data = returns)
fit2[[4]] <- felm(I(retgap*100) ~ log(tna) + I(exp_ratio*100) + turn_ratio + logmktcap + btm + mom + I(lehman_q == 5) * post + I(lehman_q == 4) * post + I(lehman_q == 3) * post + I(lehman_q == 2) * post | date | 0 | mgmt_cd, data = returns)

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE, keep.stat = c("n", "adj.rsq")) 

### Panel C ### 
rm(list = ls())
trades <- fread("pseudo_data_lehman_shortfall.csv")

fit3 <- NULL
fit3[[1]] <- felm(I(es_ticket*100) ~ quint_lehman_2007 * post 
                  + abs(ret*100) + abs(sprtrn*100) + I(side==1) + ret_lag + ret_lag * I(side==1) + log(avg_vol_lag) + I(exchcd==1) + I(1/abs(prc)) + I(volume_ticket/avg_vol_lag) 
                  | date_week + brokercode | 0 | brokercode + date_week, data = trades)
fit3[[2]] <- felm(I(es_ticket*100) ~ quint_lehman_2007 * post 
                  + abs(ret*100) + abs(sprtrn*100) + I(side==1) + ret_lag + ret_lag * I(side==1) + log(avg_vol_lag) + I(exchcd==1) + I(1/abs(prc)) + I(volume_ticket/avg_vol_lag) 
                  | date_week + brokercode + permno | 0 | brokercode + date_week, data = trades)
fit3[[3]] <- felm(I(es_ticket*100) ~ I(quint_lehman_2007 == 5) * post + I(quint_lehman_2007 == 4) * post + I(quint_lehman_2007 == 3) * post + I(quint_lehman_2007 == 2) * post  
                  + abs(ret*100) + abs(sprtrn*100) + I(side==1) + ret_lag + ret_lag * I(side==1) + log(avg_vol_lag) + I(exchcd==1) + I(1/abs(prc)) + I(volume_ticket/avg_vol_lag) 
                  | date_week + brokercode | 0 | brokercode + date_week, data = trades)
fit3[[4]] <- felm(I(es_ticket*100) ~ I(quint_lehman_2007 == 5) * post + I(quint_lehman_2007 == 4) * post + I(quint_lehman_2007 == 3) * post + I(quint_lehman_2007 == 2) * post  
                  + abs(ret*100) + abs(sprtrn*100) + I(side==1) + ret_lag + ret_lag * I(side==1) + log(avg_vol_lag) + I(exchcd==1) + I(1/abs(prc)) + I(volume_ticket/avg_vol_lag) 
                  | date_week + brokercode + permno | 0 | brokercode + date_week, data = trades)

stargazer(fit3, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE, keep.stat = c("n", "adj.rsq"))


################
### Table 10 ###
################

### Import data ### 
rm(list = ls())
holdings <- fread("pseudo_data_density_car.csv")

### Report CAR ### 
res1 <- holdings[, .(car_avg = mean(car , na.rm = TRUE)), by = .(den_q, date)][order(den_q, date)]
res2 <- dcast(res1, date ~ den_q, value.var = "car_avg")
names(res2) <- c("date", paste0("q", 1:5))

fit1 <- NULL
fit1[[1]] <- lm(q1 ~ 1, data = res2)
fit1[[2]] <- lm(q2 ~ 1, data = res2)
fit1[[3]] <- lm(q3 ~ 1, data = res2)
fit1[[4]] <- lm(q4 ~ 1, data = res2)
fit1[[5]] <- lm(q5 ~ 1, data = res2)
fit1[[6]] <- lm((q5-q1) ~ 1, data = res2)

fit2 <- NULL
for (i in 1:6) {
  fit2[[i]] <- coeftest(fit1[[i]], vcov.=NeweyWest(fit1[[i]], lag=5, adjust=FALSE, verbose=TRUE))  
}

stargazer(fit2, digits = 2, report = "vc*t", align = TRUE, no.space = TRUE)

